      subroutine bvec(nxp, nyp, nzp, ijkcell, ncells, tsix, tsiy, tsiz, etax, &
         etay, etaz, nux, nuy, nuz, rmaj, rwall, strait, toroid, q0, bzi, psi, &
         x, y, z, bx, by, bz) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nxp 
      integer , intent(in) :: nyp 
      integer , intent(in) :: nzp 
      integer , intent(in) :: ncells 
      real(double) , intent(in) :: rmaj 
      real(double)  :: rwall 
      real(double) , intent(in) :: strait 
      real(double) , intent(in) :: toroid 
      real(double) , intent(in) :: q0 
      real(double) , intent(in) :: bzi 
      real(double)  :: psi 
      integer , intent(in) :: ijkcell(*) 
      integer , intent(in) :: nux(*) 
      integer , intent(in) :: nuy(*) 
      integer , intent(in) :: nuz(*) 
      real(double) , intent(in) :: tsix(*) 
      real(double) , intent(in) :: tsiy(*) 
      real(double) , intent(in) :: tsiz(*) 
      real(double) , intent(in) :: etax(*) 
      real(double) , intent(in) :: etay(*) 
      real(double) , intent(in) :: etaz(*) 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(in) :: z(*) 
      real(double) , intent(inout) :: bx(*) 
      real(double) , intent(inout) :: by(*) 
      real(double) , intent(inout) :: bz(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: ibar, jbar, kbar, n, ijk 
      real(double) :: pi, yv, ws, wsr, wsmr, cosphi, phi, wsq0, denom, rbmag 
!-----------------------------------------------
!
      ibar = nxp - 2 
      jbar = nyp - 2 
      kbar = nzp - 2 
      pi = acos(-1.) 
!
      do n = 1, ncells 
!
         ijk = ijkcell(n) 
!
         yv = y(ijk)*toroid 
         ws = toroid*sqrt(x(ijk)**2+yv**2) - rmaj + strait*x(ijk) 
!
!     wsr is the minor radius
!
         wsr = sqrt(ws**2 + z(ijk)**2) 
!
!     wsmr is the major radius
!
         wsmr = sqrt(x(ijk)**2+y(ijk)**2+z(ijk)**2) 
!     Ref: Bauer, Betancourt, and Garabedian
!     "A Computational Method in Plasma Physics"
!
!
!
!     calculate toroidal angle
!
         cosphi = x(ijk)/wsmr 
         cosphi = amin1(1.,cosphi) 
         cosphi = amax1(-1.,cosphi) 
         phi = acos(cosphi) 
!
         wsq0 = q0*float(ibar)*cos(2.*phi)*wsr/float(jbar) 
         denom = 2.*pi/(float(ibar)*float(kbar)) 
!
         bx(ijk) = bzi*(nuy(ijk)*(wsq0*etaz(ijk)+tsiz(ijk))-nuz(ijk)*(wsq0*etay&
            (ijk)+tsiy(ijk)))*denom*(rmaj/wsmr + strait)*wsr 
!
         by(ijk) = bzi*(nuz(ijk)*(wsq0*etax(ijk)+tsix(ijk))-nux(ijk)*(wsq0*etaz&
            (ijk)+tsiz(ijk)))*denom*(rmaj/wsmr + strait)*wsr 
!
         bz(ijk) = bzi*(nux(ijk)*(wsq0*etay(ijk)+tsiy(ijk))-nuy(ijk)*(wsq0*etax&
            (ijk)+tsix(ijk)))*denom*(rmaj/wsmr + strait)*wsr 
!
         rbmag = 1./(sqrt(bx(ijk)**2+by(ijk)**2+bz(ijk)**2)+1.E-20) 
!
         bx(ijk) = bx(ijk)*rbmag 
         by(ijk) = by(ijk)*rbmag 
         bz(ijk) = bz(ijk)*rbmag 
!
      end do 
!
      return  
      end subroutine bvec 


      subroutine fluxsurf(i1, i2, j1, j2, k1, k2, iwid, jwid, kwid, x, y, z, &
         rwall, rmaj, toroid, strait, delta, psi) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: i1 
      integer , intent(in) :: i2 
      integer , intent(in) :: j1 
      integer , intent(in) :: j2 
      integer , intent(in) :: k1 
      integer , intent(in) :: k2 
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      real(double) , intent(in) :: rwall 
      real(double) , intent(in) :: rmaj 
      real(double) , intent(in) :: toroid 
      real(double) , intent(in) :: strait 
      real(double) , intent(in) :: delta 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(in) :: z(*) 
      real(double) , intent(out) :: psi(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: k, j, i, ijk 
      real(double) :: ws, rminor, xv, yv, rmr, sphi, cphi, stheta, ctheta, &
         deltar, c2theta, c3theta 
!-----------------------------------------------
!
      do k = k1, k2 
         do j = j1, j2 
            do i = i1, i2 
!
               ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid 
!
!     calculate toroidal coordinates from cartesian coordinates
!
               ws = toroid*sqrt(x(ijk)**2+y(ijk)**2) - rmaj + strait*x(ijk) 
!
               rminor = sqrt(ws**2 + z(ijk)**2) 
!
               xv = x(ijk)*toroid + strait 
               yv = y(ijk)*toroid 
!
               rmr = 1./sqrt(xv**2 + yv**2) 
!
               sphi = yv*rmr 
               cphi = xv*rmr 
!
               stheta = ws/(rminor + 1.E-35) 
               ctheta = z(ijk)/(rminor + 1.E-35) 
!
!      deltar=delta*rminor*(rwall-rminor)/rwall**2
               deltar = delta*rminor**3*(rwall - rminor)/rwall**4*256./27. 
!
!     multiple angle formulae
!
!     s2theta=2.*stheta*ctheta
               c2theta = 2.*ctheta**2 - 1. 
!
!     s3theta=3.*stheta-4.*stheta**2
               c3theta = 4.*ctheta**3 - 3.*ctheta 
!
               psi(ijk) = rminor*(1. - deltar*c2theta) 
!
            end do 
         end do 
      end do 
!
      return  
      end subroutine fluxsurf 


      subroutine refvec(i1, i2, j1, j2, k1, k2, iwid, jwid, kwid, tsix, tsiy, &
         tsiz, etax, etay, etaz, nux, nuy, nuz, x, y, z, toroid, strait, rwall&
         , rmaj, psi, cgx, cgy, cgz, agx, agy, agz, bgx, bgy, bgz) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: i1 
      integer , intent(in) :: i2 
      integer , intent(in) :: j1 
      integer , intent(in) :: j2 
      integer , intent(in) :: k1 
      integer , intent(in) :: k2 
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      real(double) , intent(in) :: toroid 
      real(double) , intent(in) :: strait 
      real(double) , intent(in) :: rwall 
      real(double) , intent(in) :: rmaj 
      real(double) , intent(in) :: tsix(*) 
      real(double) , intent(in) :: tsiy(*) 
      real(double) , intent(in) :: tsiz(*) 
      real(double) , intent(in) :: etax(*) 
      real(double) , intent(in) :: etay(*) 
      real(double) , intent(in) :: etaz(*) 
      real(double) , intent(in) :: nux(*) 
      real(double) , intent(in) :: nuy(*) 
      real(double) , intent(in) :: nuz(*) 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(in) :: z(*) 
      real(double) , intent(in) :: psi(*) 
      real(double) , intent(inout) :: cgx(*) 
      real(double) , intent(inout) :: cgy(*) 
      real(double) , intent(inout) :: cgz(*) 
      real(double) , intent(inout) :: agx(*) 
      real(double) , intent(inout) :: agy(*) 
      real(double) , intent(inout) :: agz(*) 
      real(double) , intent(inout) :: bgx(*) 
      real(double) , intent(inout) :: bgy(*) 
      real(double) , intent(inout) :: bgz(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: k, j, i, ijk, ijka 
      real(double) :: psi1, psi2, psi3, rcg, rag, ws, xv, yv, rmr, sphi, cphi, &
         stheta, ctheta 
!-----------------------------------------------
!
!
      do k = k1, k2 - 1 
         do j = j1, j2 
            do i = i1, i2 
!
               ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid 
!
               psi1 = 0.5*(psi(ijk+iwid)-psi(ijk-iwid)) 
               psi2 = 0.5*(psi(ijk+jwid)-psi(ijk-jwid)) 
               psi3 = 0.5*(psi(ijk+kwid)-psi(ijk-kwid)) 
!
               cgx(ijk) = psi1*tsix(ijk) + psi2*etax(ijk) + psi3*nux(ijk) 
               cgy(ijk) = psi1*tsiy(ijk) + psi2*etay(ijk) + psi3*nuy(ijk) 
               cgz(ijk) = psi1*tsiz(ijk) + psi2*etaz(ijk) + psi3*nuz(ijk) 
!
               rcg = 1./(sqrt(cgx(ijk)**2+cgy(ijk)**2+cgz(ijk)**2)+1.E-30) 
!c
               cgx(ijk) = cgx(ijk)*rcg 
               cgy(ijk) = cgy(ijk)*rcg 
               cgz(ijk) = cgz(ijk)*rcg 
!
!     construct a vector that is orthogonal to b (calculated in bvec)
!     and c, which is normal to flux surfaces
!
               agx(ijk) = bgy(ijk)*cgz(ijk) - bgz(ijk)*cgy(ijk) 
               agy(ijk) = bgz(ijk)*cgx(ijk) - bgx(ijk)*cgz(ijk) 
               agz(ijk) = bgx(ijk)*cgy(ijk) - bgy(ijk)*cgx(ijk) 
!
               rag = 1./(sqrt(agx(ijk)**2+agy(ijk)**2+agz(ijk)**2)+1.E-20) 
!
               agx(ijk) = agx(ijk)*rag 
               agy(ijk) = agy(ijk)*rag 
               agz(ijk) = agz(ijk)*rag 
!
            end do 
         end do 
      end do 
!
!     calculate reference vectors on the axis and the outer wall
!
      do j = j1, j2 
         do i = i1, i2 
!
            ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (k2 - 1)*kwid 
!
            ws = toroid*sqrt(x(ijk)**2+y(ijk)**2) - rmaj + strait*x(ijk) 
!
            xv = x(ijk)*toroid + strait 
            yv = y(ijk)*toroid 
!
            rmr = 1./sqrt(xv**2 + yv**2) 
!
            sphi = yv*rmr 
            cphi = xv*rmr 
            stheta = ws/rwall 
            ctheta = z(ijk)/rwall 
!
            cgx(ijk) = sphi*ctheta 
            cgy(ijk) = cphi*ctheta 
            cgz(ijk) = stheta 
            agx(ijk) = agx(ijk-kwid) 
            agy(ijk) = agy(ijk-kwid) 
            agz(ijk) = agz(ijk-kwid) 
!
            bgx(ijk) = bgx(ijk-kwid) 
            bgy(ijk) = bgy(ijk-kwid) 
            bgz(ijk) = bgz(ijk-kwid) 
!
!      cgx(ijk)=cgx(ijk-kwid)
!      cgy(ijk)=cgy(ijk-kwid)
!      cgz(ijk)=cgz(ijk-kwid)
!
!
!     set unit vectors on the axis equal to unit vectors at the wall
!
            ijka = 1 + (i - 1)*iwid + (j - 1)*jwid 
!
            cgx(ijka) = cgx(ijk) 
            cgy(ijka) = cgy(ijk) 
            cgz(ijka) = cgz(ijk) 
!
!
         end do 
      end do 
!
!
      return  
      end subroutine refvec 
